11/30/2018
Probability model.
Model encodes our understanding of the scientific process of interest.
Model accounts for as much uncertainty as possible.
Model results in a probability distribution.
Update model with data.
Criticize the model
Does the model fit the data well?
Do the predictions make sense?
Are there subsets of the data that don't fit the model well?
Make inference using the model.
Start with probability distributions:
\[ \begin{align*} \left[y_i | \boldsymbol{\theta} \right] & \sim \operatorname{N}(X_i \beta, \sigma^2) \\ \boldsymbol{\theta} & = (\beta, \sigma^2) \end{align*} \]
Hierarchical model:
A model built in components.
Each component represents a different statistical goal.
Break the model into components:
Data Model.
Process Model.
Prior Model.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} [\mathbf{z} | \boldsymbol{\theta}_P] [\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} }% \]
\[ {\huge \begin{align*} \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} \end{align*} } \]
\[ {\huge \begin{align*} \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} \end{align*} } \]
\[ {\huge \begin{align*} \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} \end{align*} } \]
Age of minerals:
\(\mathbf{y}\) is the radio-date estimate.
\(\mathbf{z}\) is the true mineral age.
\(\theta_D\) is the radio-date standard error.
The probability distribution is determined by the measurement process.
\[ {\huge \begin{align*} \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} \end{align*} } \]
Reconstructing climate from tree rings
\(\mathbf{y}\) is the tree ring width increment.
\(\mathbf{z}\) is the true, unobserved climate variable.
\(\boldsymbol{\theta}_D\) models the relationship between climate, stand dynamics, individual heterogeneity, tree age, (etc.) and tree ring width.
The probability distribution is determined by tree physiology, measurement uncertainty, etc.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]}[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
\[ {\huge \begin{align*} \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]} \end{align*} } \]
\[ {\huge \begin{align*} \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]} \end{align*} } \]
\[ {\huge \begin{align*} \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]} \end{align*} } \]
Reconstructing climate with tree rings
Trees of the same species share a similar response to climate.
Climate variables at sites nearby in location are closer to each other than sites far apart, on average.
Climate variables seperated by short periods of time are more similar than climate variables over long periods of time.
\(\mathbf{z}\) is the value of the unobserved climate variables.
\(\boldsymbol{\theta}_P\) are the species-specific growth responses and the correlation of climate across time and space.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] [\mathbf{z} | \boldsymbol{\theta}_P] \color{orange}{[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P]} \end{align*} } \]
Probability distributions define "reasonable" ranges for parameters.
\[ {\huge \begin{align*} \color{cyan}{[\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}]} & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] [\mathbf{z} | \boldsymbol{\theta}_P] [\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
\[ {\huge \begin{align*} \color{cyan}{[\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}]} \end{align*} } \]
Probability distribution over all unknowns in the model.
Inference is made using the posterior distribution.
Because the posterior distribution is a probability distribution, uncertainty is easy to calculate.
Climate change is well understood globally.
Climate change is less well understood locally.
Need for spatially explicit reconstructions of climate variables.
Problem: data sources are messy and noisy.
By OakleyOriginals - https://www.flickr.com/photos/oakleyoriginals/2954878643/, CC BY 2.0, Link
Vegetation composition and structure change from ice age to current period.
Using change in temperature to predict future vegetation change.
Data model: Multi-logit distribution for ordered categories of observed change.
Process model: Assumes increasing temperature results in smooth changes of composition and struction.
Prior model: Not used.
Sharman and Johnstone (2017). Sediment unmixing using detrital geochronology. Earth and Palenetary Science Letters.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} [\mathbf{z} | \boldsymbol{\theta}_P] [\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} }% \]
\[ \begin{align*} y_{ib} & \sim \color{red}{\operatorname{N}(z_{ib}, \sigma^2_{ib}).} \\ y_{id} & \sim \color{red}{\operatorname{N}(z_{id}, \sigma^2_{id}).} \end{align*} \]
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]}[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
Assumptions:
\[ \begin{align*} {z}_{ib} \sim \color{blue}{\sum_{k=1}^K p_{bk} \operatorname{N}(\mu_k, \sigma^2_k).} \end{align*} \]
\[ \begin{align*} z_{id} & \sim \color{blue}{\sum_{b=1}^B \phi_b \sum_{k=1}^K p_{bk} \operatorname{N}(\mu_k, \sigma^2_k).} \end{align*} \]
\[ \begin{align*} \phi_1 = 0.200 \quad\quad\,\,\, \phi_2 = 0.532 \quad\quad\,\,\,\,\, \phi_3 = 0.268 \,\,\quad\quad\quad \mbox{Daughter} \end{align*} \]
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] [\mathbf{z} | \boldsymbol{\theta}_P]\color{orange}{[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P]} \end{align*} } \]
\[ \begin{align*} y_{id} & \sim \color{red}{\operatorname{N}(z_{id}, \sigma^2_{id}).} \end{align*} \]
\[ \begin{align*} z_{id} & \sim \color{blue}{\sum_{b=1}^B \phi_{db} \sum_{k=1}^K p_{bk} \operatorname{N}(\mu_k, \sigma^2_k).} \end{align*} \]
\(\color{orange}{\mbox{Prior model:}}\)
\[ \begin{align*} \sum_{k=1}^K I\{ \phi_b^{(k)} > 0.5 \}. \end{align*} \]
Account for spatial correlation among daughters.
Account for temporal correlation within a sediment core.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} [\mathbf{z} | \boldsymbol{\theta}_P] [\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} }% \]
\[ \begin{align*} \begin{pmatrix} \mathbf{T}_{t} \\ \mathbf{P}_{t} \end{pmatrix} = \begin{pmatrix} \mathbf{T}_{\mbox{Jan}, t} \\ \vdots \\ \mathbf{T}_{\mbox{Dec}, t} \\ \mathbf{P}_{\mbox{Jan}, t} \\ \vdots \\ \mathbf{P}_{\mbox{Dec}, t} \\ \end{pmatrix} \sim \color{red}{\mbox{N} \left( \mathbf{A} \begin{pmatrix} \mathbf{T}_{t-1} \\ \mathbf{P}_{t-1} \end{pmatrix}, \boldsymbol{\Sigma} \right).} \end{align*} \]
\[ \begin{align*} y_{i t j} & \sim \color{red}{\begin{cases} \mbox{N}\left(\beta_{0_j} + \beta_{1_j} f^{VS}\left(\mathbf{w}_t, \boldsymbol{\theta}^{VS}_j\right), \sigma^2_j \right) & \mbox{if } z_j = 0,\\ \mbox{N}\left(\tilde{\beta}_{0_j} + \tilde{\beta}_{1_j} f^{Pro}\left(\mathbf{w}_t, \boldsymbol{\theta}^{Pro}_j\right), \tilde{\sigma}^2_j \right) & \mbox{if } z_j = 1. \\ \end{cases}} \end{align*} \]
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]}[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
\[ \begin{align*} f^{\ell}\left(\mathbf{w}_t, \boldsymbol{\theta}^{\ell}_j\right) & = \color{blue}{ \sum_{s=1}^{12} \chi_s \mbox{min} \left( g^{\ell}\left( T_{t,s}, \boldsymbol{\theta}^{\ell}_j \right), g^{\ell}\left( P_{t, s}, \boldsymbol{\theta}^{\ell}_j \right) \right),} \\ & \ell = VS \mbox{ or } Pro. \end{align*} \]
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] [\mathbf{z} | \boldsymbol{\theta}_P]\color{orange}{[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P]} \end{align*} } \]
Assumes growth responses follow "ecological niche."
Tree species that grow in the Hudson Valley respond to similar climate so have similar responses.
Variations from common response are to exploit an "ecological niche" that allows many species to exist on the same landscape.
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto \color{red}{[\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}]} [\mathbf{z} | \boldsymbol{\theta}_P] [\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
\[ \begin{align*} \mathbf{y}\left( \mathbf{s}_i, t \right) & \sim \color{red}{\operatorname{Dirichlet-Multinomial} \left( N, \exp \left(\mathbf{z}\left( \mathbf{s}_i, t \right) \boldsymbol{\beta} \right)\right)} \end{align*} \]
Vegetation response to climate is non-linear.
Pollen are "aggregated" into groups across space and taxonomy.
\[ \begin{align*} \mathbf{y}\left( \mathbf{s}_i, t \right) & \sim \color{red}{\operatorname{Dirichlet-Multinomial} \left( N, \exp\left( \mathbf{B} \left( \mathbf{z}\left( \mathbf{s}_i, t \right) \right) \boldsymbol{\beta} \right) \right)} \end{align*} \]
\[ {\huge \begin{align*} [\mathbf{z}, \boldsymbol{\theta}_D, \boldsymbol{\theta}_P | \mathbf{y}] & \propto [\mathbf{y} | \boldsymbol{\theta}_D, \mathbf{z}] \color{blue}{[\mathbf{z} | \boldsymbol{\theta}_P]}[\boldsymbol{\theta}_D] [\boldsymbol{\theta}_P] \end{align*} } \]
\[\ \begin{align*} \color{blue}{\mathbf{z} \left(t \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma}} & \sim \color{blue}{\operatorname{N}\left(\mathbf{M}\left(t\right) \left( \mathbf{z} \left(t-1 \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} \right), \boldsymbol{R}\left( \boldsymbol{\theta} \right) \right)} \end{align*} \]
Banerjee S, Gelfand AE, Finley AO and Sang H (2008). "Gaussian, predictive process models for large spatial data sets." Journal, of the Royal Statistical Society: Series B (Statistical, Methodology), 70(4), pp. 825-848.
\(\boldsymbol{\eta} \left( t \right) \approx \mathbf{r} \left(\mathbf{s}, \mathbf{s}^{\star} \right) \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)^{-1} \tilde{\boldsymbol{\eta}} \left( t \right)\).
The predictive process can be shown to be the optimal predictor of the parent process \(\boldsymbol{\eta} \left( t \right)\) of dimension \(m\)
The dynamic climate process becomes
\(\begin{align*} \mathbf{z} \left(t \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} & \approx \mathbf{M}\left(t\right) \left( \mathbf{z} \left(t-1 \right) - \mathbf{X} \left( t \right) \boldsymbol{\gamma} \right) + \mathbf{r} \left(\mathbf{s}, \mathbf{s}^{\star} \right) \mathbf{R}^{\star}\left( \boldsymbol{\theta} \right)^{-1} \boldsymbol{\tilde{\eta}} \left(t \right) \end{align*}\)
devtools::install_github("jtipton25/BayesComposition")
C++ using Rcpp package for computation speed.Model framework opens the door to answering meaningful questions:
It is possible to put the science in your statistics.
Takes some careful thinking and learning.
Opens the door to more powerful analyses.
More flexibility in the questions that can be answered.